Microconfined High-Pressure Transcritical Channel Flow Database: Laminar, Transitional & Turbulent Regimes

The potential of comprehending and managing microscale flows to enhance energy processes, especially in heat transfer and propulsion applications, remains largely untapped particularly for supercritical fluids, which have gained increased interest over the past years due to the higher power and thermodynamic efficiencies they provide. This work, therefore, presents the first comprehensive, open-source dataset carefully curated and structured for studying microconfined high-pressure transcritical fluid channel flows under various regimes. Particularly, the dataset contains 18 direct numerical simulations of carbon dioxide at different bulk pressures and velocities confined between differentially-heated walls. For all cases, the thermodynamic conditions selected impose the fluid to undergo a transcritical trajectory across the pseudo-boiling region. The data collection comprises an array of physical quantities that enable comprehensive parametric analyses spanning laminar, transitional, and turbulent flow regimes. This data repository is poised to provide access to the detailed study and modeling of the complex flow physics observed in high-pressure transcritical fluids, especially those closely linked to improving microfluidics performance.


Background & Summary
Over the past decades, the field of microfluidics has experienced rapid expansion driven by its extraordinary high surface-to-volume ratios, precise flow control, and well-matched length scales for interacting with microscopic elements 1 .However, typical microfluidic systems face limitations in terms of flow mixing and heat transfer due to their inherent operation under laminar flow regimes as a result of their characteristic small dimensions (H ∼1−1000 μm) and relatively low bulk velocities ( < U 1 b m/s).Specifically, at atmospheric pressure conditions, the bulk Reynolds numbers b encountered in microfluidics fall within the approximate range 0.1-100, where the bulk kinematic viscosity b ν typically ranges from 10 −6 to 10 −4 m 2 /s.In this regard, mixing within microdevices stands as one of the pivotal technologies for enabling miniaturized analysis and mass & energy transfer/conversion solutions.For instance, selected microfluidic applications in which their performance depends upon how rapidly homogeneous mixing is generated include high-throughput chemical synthesis, preparation of emulsions, protein folding, flash chemistry, and high-pressure chromatography.In this context, a crucial field of study in microfluidics revolves around enhancing molecular diffusion within laminar flows through various strategies, like the fabrication of microchannels with serpentine patterns 2 .Nonetheless, this technique frequently introduces significant complexities in manufacturing, without delivering remarkable increases in mixing and transfer rates.Alternative strategies that have been extensively investigated, for example in the works of Wang et al. 3 & Nan et al. 4 , attempt to enhance chaotic mixing through electrokinetic forcing.This type of approaches, however, requires the utilization of fluids with significantly different electric properties.A different (brute-force) strategy is to extraordinarily increase the volumetric flow rate (large sizes/velocities) in microdevices to reach incipient turbulent flow conditions [5][6][7] .Although interesting for some applications, for instance high throughput liquid-liquid extraction 8 , this approach presents two main drawbacks: (i) there is no intrinsic mechanism to trigger flow destabilization, and (ii) the volumetric flow rates required are notably high in contradiction to the small rates typically sought in microfluidic applications.Finally, an interesting approach recently proposed is the utilization of cross-flow micro-jets 9 to enhance the generation of large-scale (laminar) vortices in microfluidic systems, which aids in increasing mixing and heat transfer.
Nonetheless, a different promising strategy, which differs from conventional microfluidic practices, has been recently proposed 10,11 based on inducing microconfined turbulence by operating under high-pressure transcritical conditions to harness the unique thermophysical properties of supercritical fluids.Supercritical fluids, which exist at pressures above the critical pressure P c , undergo a complex second-order phase transition charac- terized by substantial property changes across the pseudo-boiling region, i.e., transcritical trajectory 10,12 .This transition, which occurs at the pseudo-boiling temperature T pb and is responsible for flow destabilization through baroclinic torque effects 13 , divides the fluid into two distinct regions marked by non-ideal behavior.Specifically, this behavior is characterized by significant deviations from the ideal-gas law, where the compressibility factor Z deviates from unity, caused by the modification of the molecular structure of the fluid resulting from the phase change.In detail, at temperatures below the pseudo-boiling temperature ( < T T pb ), the fluid behaves as a supercritical liquid-like state characterized by: (i) high density ρ; (ii) low isothermal compressibility T T ; and (iii) decreasing dynamic viscosity μ and thermal conductivity κ with temperature.Conversely, when the temperature exceeds the pseudo-boiling point ( > T T pb ), the fluid is at supercritical gas-like state and exhibits opposite characteristics/behaviors. Furthermore, as discussed in previous studies 12,14 , it is worth noting that the utilization of high-pressure conditions is not limited to microfluidic research but also extends to a broad spectrum of engineering applications, such as supercritical water-cooled reactors and liquid rocket engines.This underscores the broader significance of the unconventional method discussed earlier and its potential impact on diverse technological solutions.In this regard, studying supercritical fluids experimentally poses considerable challenges, primarily because of the elevated pressures typically required, which complicate optical and sensing access.The thermophysical complexity of such systems is further heightened by the interest to operate under microfluidic conditions.Hence, the accurate characterization and prediction of high-pressure transcritical turbulent flows typically necessitate the use of computational scale-resolving approaches, such as direct numerical simulation (DNS) and large-eddy simulation (LES).Several notable examples of scale-resolving studies in the field include DNS investigations of supercritical heat transfer in pipe flows conducted by Bae et al. 15 , studies of transcritical flows near the critical point carried out by Sengupta et al. 16 , and research by Kawai 17 involving DNS of transcritical turbulent boundary layers.
In the context of microconfined supercritical fluids, Bernades et al. 13 conducted a DNS of channel flow using Nitrogen as the working fluid.Differently to the thermodynamic regimes considered in this work, the study focused only on 2 specific differentially-heated scenarios at atmospheric (ideal gas) and high-pressure (real gas) conditions.The simulations aimed to achieve a friction Reynolds number of Re 100 = τ at the cold/bottom wall to explore the possibility of inducing turbulence under these conditions.In this regard, motivated by the lack of sufficient research data and the escalating scientific interest on high-pressure transcritical flows at microfluidic conditions, the present work aims to carefully compute a set of DNSs to generate a curated database of 18 microconfined high-pressure transcritical channel cases encompassing laminar, transitional, and turbulent flow regimes.To that end, the numerical methods utilized for the simulations are presented, along with a comprehensive description of the framework employed for the modeling of supercritical fluids.In doing so, the first curated and logically structured computational microconfined high-pressure transcritical channel flow dataset is presented for immediate use, enabling researchers to analyze and explore the intricacies of such complex flows, along with the availability of data for the replication of computational experiments.

Methods
The mathematical modeling utilized for studying microconfined supercritical fluids in terms of (i) equations of fluid motion, (ii) real-gas thermodynamics, (iii) high-pressure transport coefficients, and (iv) numerical methods is described below.This framework has been previously utilized to investigate the physics of high-pressure transcritical flows 11,13,18,19 , develop advanced numerical schemes 20,21 , and perform data-driven analysis 22,23 .The section culminates with an in-depth exposition of the specific computational experiments featured in the dataset.
equations of fluid motion.The motion of high-pressure transcritical fluids is described by the conservation of mass, momentum, and total energy, which in dimensionless form are written as where superscript denotes normalized quantities, t is the time, u is the velocity vector, ρ is the density, P is the pressure, ) is the viscous stress tensor with μ the dynamic viscosity and I the identity matrix, E e u /2 2 = + and e are the total and internal energy, respectively, T is the temperature, and κ is the thermal conductivity.The obtention of these dimensionless equations is fundamented on the following set of inertial-based scalings 23,24  accounting for the ratio between advective mass transfer and heat dissipation potential.
Real-gas thermodynamics.The thermodynamic space of solutions for the state variables pressure P, temperature T, and density ρ of a monocomponent fluid is described by an equation of state.In this work, therefore, the Peng-Robinson equation of state 25 , which is suitable for high-pressure thermodynamic conditions as provided in the work by Jofre & Urzay 12 , is utilized to close the equations of fluid motion above.In general form, it can be expressed in terms of the compressibility factor the specific gas constant of CO 2 , which in dimensionless form is written as is an approximated real-gas heat capacity ratio 26 with c V the isochoric heat capacity.As it can be noted, the dimensionless bulk Mach number Ma U c / b b b = appears, where c b is the bulk speed of sound, which represents the ratio of flow velocity to the local speed of sound.In addition, real-gas equations of state need to be supplemented with the corresponding high-pressure thermodynamic variables (e.g., internal energy, heat capacities) based on departure functions 27 calculated as a difference between two states.In particular, their usefulness is to transform thermodynamic variables from ideal-gas conditions (low pressure -only temperature dependant) to supercritical conditions (high pressure).The ideal-gas parts are calculated by means of the NASA 7-coefficient polynomial 28 , while the analytical departure expressions to high pressures are derived from the Peng-Robinson equation of state 12 .
High-pressure transport coefficients.The high pressures involved in the analyses conducted in this work prevent the use of simple relations for the calculation of dynamic viscosity μ and thermal conductivity κ.In this regard, standard methods for computing these coefficients for Newtonian fluids are based on the correlation expressions proposed by Chung et al. 29,30 .These correlation expressions are mainly function of critical temperature T c and density ρ c , molecular weight W, acentric factor ω, association factor κ a and dipole moment M, and the NASA 7-coefficient polynomial 28 ; further details can be found in dedicated works, like for example Poling et al. 31 and Jofre and Urzay 12 .
Numerical method.The equations of fluid motion are computationally solved by means of the in-house flow solver RHEA 32 .A standard semi-discretization procedure is adopted in which they are firstly discretized in space and then integrated in time.In particular, spatial operators are treated using second-order central-differencing schemes, and time-advancement is explicitly performed by means of a third-order strong-stability preserving (SSP) Runge-Kutta approach 33 coupled with an artificial compressibility method (ACM) 34 .The convective terms are expanded according to the Kennedy-Gruber-Pirozzoli (KGP) splitting 35,36 , which has been recently extended to high-pressure transcritical fluids turbulence in previous works 20,37,38 .The method preserves kinetic energy by convection, and is locally conservative for mass, momentum, and total energy.This numerical framework provides stable computations without the need of any form of artificial dissipation or stabilization procedures.

Computational experiments.
As portrayed in Fig. 1, a high-pressure transcritical microfluidic channel flow at low-Mach-number conditions is chosen for the dataset generation.The selected fluid is CO 2 with its critical thermophysical conditions and properties outlined in Table 1.The microchannel is operated at three distinctive supercritical bulk pressures = .
, 2 and 5, with the fluid enclosed between two parallel isothermal walls at cold (cw) and hot (hw) temperatures expressed as T cw /T c and T hw /T c , respectively.These temperature values lie within the specified intervals ΔT/T c = (T hw − T cw )/T c = 0.8-1.4,0.9-1.2 and 0.95-1.1.The separation between walls is equal to H = 2δ, with δ = 100 μm representing half of the channel height.The fluid flows from left to right in the streamwise direction with bulk velocities U b = 0.5 and 1 m/s, which is imposed through a body force controlled by a proportional feedback loop with gain k p = 0.1 aimed at reducing the difference between the desired and measured (numerical) values; the relative errors achieved are below 0.5% for all cases.This range of parameters compels the fluid to undergo varying flow conditions as it traverses the pseudo-boiling line through a transcritical thermodynamic trajectory 13,22 .In particular, as detailed in Table 2 and qualitatively represented in Fig. 2, the resulting dataset comprises a total of 18 simulations encompassing laminar (non-symmetric), transitional and turbulent flow regimes 39 .It is important to highlight that, although it is difficult to observe from the streamwise velocity colormaps of Fig. 2, none of the cases presented are trivially symmetric as they present spanwise oscillations/fluctuations of flow quantities and fluid properties.In detail, the baroclinic-type instabilities generated at the pseudo-boiling region modulate the flow in the three directions, and consequently no complete symmetry is observed for any case.
The computational domain is 4 2 (4/3) πδ δ π δ × × in the streamwise (x), wall-normal (y), and spanwise (z) directions, respectively, which is large enough to represent the largest flow scales of the problem 40 .The streamwise and spanwise boundaries are set periodic, and no-slip conditions are imposed on the horizontal boundaries (x-z planes).The resulting DNS grid contains 96 × 96 × 96 points uniformly distributed in the streamwise and spanwise directions.In particular, the mesh resolutions in cw-units fall within the ranges .< Δ < .+ from the hot wall, with the latter being the most critical one; extensive details of the resolutions in wall units can be found in the technical validation section.For the temporal time stepping, a Courant-Friedrichs-Lewy number of CFL = 0.1 is selected.For each case, the computational strategy starts from a linear velocity profile with random fluctuations 41

Data Records
The dataset is available online on Figshare 40 , a commonly used platform for sharing large datasets.The data is structured in a manner where individual cases are allocated in dedicated folders, facilitating in this way convenient access and data processing.Each case folder is structured in different FTT subfolders, which contain 10 snapshot files of the corresponding FTT uniformly spaced in time.Additionally, there is a README.txt file in each case folder that provides essential information about the individual case, such as (i) key independent parameters, (ii) dimensionless numbers, (iii) details on the calculation of the stretching factor in the y-direction, and (iv) the time advancement technique.The dataset files within the FTT subfolders are stored in the Hierarchical Data Format version 5 (h5) 42 , an open-source file format capable of handling complex and diverse data.There is also an associated xdmf file that contains data descriptions for these files, which can be visualized with, for example, the open-source software Paraview 43 .Finally, each case folder also includes a configuration file and a substances library file in YAML format.Likewise, there are C++ files in both .cppand .hppformats, allowing users to replicate the simulations for further investigation and/or verification.The dataset includes an extensive array of input features and corresponding labels, as outlined in Table 2.It has been meticulously organized to facilitate the extraction of two distinct categories of data: (i) base output variables, and (ii) derived quantities.The base output variables entail critical parameters that are directly generated by the simulations, serving as essential prerequisites for the computation of other derived quantities.These base outputs primarily consist of instantaneous flow parameters.In contrast, the derived quantities comprise various statistical measures, including first-and second-order Reynolds-and Favre-averaged statistical Tt dt ( ) Derived time-averaged quantities of the database cases.
quantities.For a given quantity, denoted as φ, the Reynolds-averaged mean φ and its fluctuation φ ′ are defined as φ φ φ = + ′ ¯, with the condition that the average of φ ′ equals zero.Similarly, the Favre-averaged mean φ and its fluctuation φ ′′ are defined as φ φ φ = + ′′ with φ ρφ ρ = ¯ / .A summary of the primary base output fields available within the dataset is provided in Table 3, while Tables 4, 5, and 6 furnish concise overviews of the derived output fields.

Technical Validation
The computational framework utilized in this work is comprehensively described by Bernades et al. 20 .Additionally, the flow solver RHEA 32 has been extensively validated against canonical flows at ideal-gas and high-pressure transcritical conditions in the works by Abdellatif, Barea, Masclans & Monteiro et al. 18,19,22,34 .In particular, the high-pressure thermophysical framework implemented in RHEA 32 , comprising the Peng-Robinson 25 real-gas equation of state and thermodynamic potentials, as well as the transport coefficients modeled by Chung et al. 29,30 , is validated through a comparison with the intermediate pressure scenario of the dataset cases at P P / 2 c = .This validation is conducted across a range of temperatures that are available in the , closely matches the reference data for density and dynamic viscosity.Notably, the thermophysical model effectively captures the rapid variations in these properties across the pseudo-boiling region at T T / 11 pb c ≈ . ,which corresponds to a second-order phase transition from supercritical liquid-like to gas-like fluid.
Additionally, when creating the dataset, a key consideration is achieving adequate mesh resolution in the three spatial directions.In this regard, as shown in Fig. 4, the grid resolution is comparable to the Kolmogorov scale η u , which is determined as η .It is important to note that the resolution values fall significantly below the accepted ranges for channel flow with highly temperature-dependent thermophysical properties [45][46][47] , which suggests that resolutions of η η Δ < are sufficient for capturing all the relevant flow scales.In particular, for the low Reynolds numbers considered in this database, by utilizing resolutions that are finer than the ones employed by reference research works, like for example Yang et al. 48, it is ensured that even high-intensity wall-shear stress events are resolved with 99% of accuracy.Likewise, Kim et al. 49 mentioned that Δ = .

Usage Notes
A post-processing script designed to facilitate the output of main variables alongside the computation of derived quantities is available in the main folder.The diligently curated dataset is perfectly tailored for direct use, providing easily accessible quantities extracted from the HDF5 files of each snapshot.Each field within the dataset constitutes a 3D array comprising 98 × 98 × 98 data points.It is worth noting that the initial and final cells in every direction, designated as boundary cells, are reserved for wall calculations.Additionally, users possess the flexibility to navigate the mesh in accordance with their individual computational needs, with the understanding that the mesh coordinates are denoted as (k, j, i), where the x, y, and z coordinates correspond to the i, j, and k directions respectively, as depicted in Fig. 5.An illustrative example of harnessing the dataset to formulate a straightforward computation of the bulk Reynolds number Re b can be found in Fig. 6.This example serves as a practical demonstration of how to effectively employ various field formats and navigate through the multiple directions of the mesh.In this illustrative case, it is imperative to note that, owing to the implementation of a stretching discretization methodology at non-slip boundaries, the precise computation of the bulk value requires the application of a volume-weighted average to the relevant quantities.
b indicating bulk quantities, x the position vector, H the total channel height, and U the time-averaged streamwise velocity.The resulting set of scaled equations includes three dimensionless numbers: c P is the isobaric heat capacity, quantifying the ratio between momentum and thermal diffusivity; and (iii) bulk Eckert number = while the mesh is stretched towards the walls in the wall-normal direction such that the first grid point is at y yu , which is advanced in time for approximately 6 and 11 flow-through-time (FTT) units for the U b = 0.5 m/s and U b = 1 m/s cases, respectively, after achieving statistically steady-state flow conditions (approximately 10 FTTs); based on the bulk velocity U b and the length of the channel L 4x πδ = , a FTT is defined as t

Fig. 1
Fig. 1 Schematic of a microconfined channel flow domain with cold (cw) and hot (hw) isothermal walls.The illustration also shows the domain dimensions and the mean flow direction.

Fig. 2
Fig. 2 Snapshots of instantaneous streamwise velocity in wall units μ + on a z/δ-δ y/ slice for the 18 cases.

Fig. 3 Fig. 4 Fig. 5
Fig. 3 Comparison between NIST data 44 and RHEA results of density normalized by critical density ρ ρ / c (a) and dynamic viscosity normalized by critical dynamic viscosity μ μ / c (b) for CO 2 at = P P / 2 c in the range ≈ .− .T T / 08 1 6 c

Fig. 6
Fig. 6 Snippet from the Python code accessing data and calculating the bulk Reynolds number from output data files.
cient to ensure third-order convergence of flow statistics.Similarly, Li et al.50 utilized a mesh similar to the one in this study to investigate up to fourth-order moments of fluctuations for high-pressure transcritical fluid flows with

Table 1 .
44itical pressure, temperature and density, molecular weight and acentric factor for CO 2 from NIST data44.

Table 2 .
Input parameters and dimensionless numbers of the database cases.

Table 3 .
Base output quantities of the database cases.

Table 5 .
Derived Reynolds-averaged root-mean-square fluctuations of the database cases.

Table 6 .
Derived Favre-averaged fluctuations of the database cases.